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ABSTRACT 

Context. 

Aims. The cosmic microwave background (CMB) power spectrum is a powerful cosmological probe as it entails almost 
all the statistical information of the CMB perturbations. Having access to only one sky, the CMB power spectrum 
measured by our experiments is only a realization of the true underlying angular power spectrum. In this paper we aim 
to recover the true underlying CMB power spectrum from the one realization that we have without a need to know the 
cosmological parameters. 

Methods. The sparsity of the CMB power spectrum is first investigated in two dictionaries; Discrete Cosine Transform 
(DCT) and Wavelet Transform (WT). The CMB power spectrum can be recovered with only a few percentage of the 
coefficients in both of these dictionaries and hence is very compressible in these dictionaries. 

Results. We study the performance of these dictionaries in smoothing a set of simulated power spectra. Based on this, 
we develop a technique that estimates the true underlying CMB power spectrum from data, i.e. without a need to know 
the cosmological parameters. 

Conclusions. This smooth estimated spectrum can be used to simulate CMB maps with similar properties to the true 
CMB simulations with the correct cosmological parameters. This allows us to make Monte Carlo simulations in a given 
project, without having to know the cosmological parameters. The developed IDL code, TOUSI, for Theoretical pOwer 
spectrUm using Sparse estimation, will be released with the next version of ISAP. 
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1. Introduction 

Measurements of the CMB anisotropics are powerful cosmological probes. In the currently favored cosmological model, 
with the nearly Gaussian-distributed curvature perturbations, almost all the statistical information are contained in the 
CMB angular power spectrum. The observed quantity on the sky is generally the CMB temperature anisotropy Q(p) in 
direction p, which is described as T(p) = TcmbW + @(p)]- This field is expanded on the spherical harmonic functions as 

+oo I 

6 (p)=E E a[l,m]Yim(p) , (1) 

1=0 m=-t 

where a[£,m] = f Q(p)Yl m (f)dp , (2) 

S 2 C M 3 is the unit sphere, £ is the multipole moment which is related to the angular size on the sky as t ~ 180° /6 and 
m is the phase ranging from —I to £. The a[£, m] are the spherical harmonic coefficients of the (noise-free) observed sky. 
For a Gaussian random field, the mean and covariance are sufficient statistics, meaning that they carry all the statistical 
information of the field. In case where the random field has zero mean, E(aoo) = and the expansion can be started 
at £ = 2, neglecting the dipole terms, i.e. £ = l 1 . For £ 2, the triangular array (a[£,m])e, m represents zero-mean, 
complex-valued random coefficients, with variance 

E(\a[£,m}\ 2 ) = C[£] > , (3) 
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1 The dipole anisotropy is dominated by the Earth's motion in space and it is hence ignored. 



2 



P. Paykari et al.: True CMB Power Spectrum Estimation 



where C[£] is the CMB angular power spectrum, which only depends on £ due the isotropy assumption. Therefore, from 
(3), an unbiased estimator of C[£] is given by the empirical power spectrum 

m 

Furthermore, as the random field is stationary, the spherical harmonic coefficients are uncorrelated, 

E{a[e,m]a*[i',m'}) = 5 lt 5 mm ,C[£] . (5) 

Since they are Gaussian they are also independent. The angular power spectrum depends on the cosmological parameters 
through an angular transfer function Tg(k) as 

C[£\=in J — T^{k)P(k) , (6) 

where k defines the scale and P{k) is the primordial matter power spectrum. 

Making accurate measurements of this power spectrum has been one of the main goals of cosmology in the past two 
decades. We have seen a range of ground- and balloon-based experiments, such as Acbar (Reichardt et al. 2009) and 
CBI (Rcadhead et al. 2004), as well as satellite experiments, such as WMAP (Bennett et al. 2003) and the recently 
launched satellite Planck (Planck Collaboration et al. 2011). All these experiments produce a temperature map of the 
sky from which the CMB power spectrum is obtained. The estimation of the power spectrum from CMB experiments 
is of great importance, as this spectrum is a way to estimate the cosmological parameters that describe the Universe. 
General methods for extracting this spectrum from a Ap^-map, with nonuniform coverage and correlated noise are quite 
expensive and time-consuming. Especially in the case of Planck where we will be dealing with N p i X = 5 x 10 7 . All these 
experiments estimate the CMB angular power spectrum from a sky map, which is a realization of the underlying true 
power spectrum; no matter how much the experiments improve, we are still limited to an accuracy within the cosmic 
variance. This means that even if we had a perfect experiment (i.e. with zero instrumental noise) we would not be able 
to recover a perfect power spectrum due to the cosmic variance limit. 

In this paper we investigate the possibility of estimating the true underlying power spectrum from a realized spectrum; 
an estimation of the true power spectrum without a need to know the cosmological parameters. For this we exploit the 
sparsity properties of the CMB power spectrum, and capitalize on it to propose an estimator of the theoretical power 
spectrum. 

The idea of sparsity in different dictionaries has previously been used. For example, (Mukhcrjec & Wang 2004) use 
wavelets to estimate the level of non-Gaussianity in the first year WMAP data. In (Mukherjee & Wang 2003) wavelets 
were used to estimate the primordial power spectrum. (Fay et al. 2008) have used wavelets to estimate the power 
spectrum from a CMB map, as an alternative to the MASTER method of (Hivon et al. 2002). There has also been 
previous attempts to smooth the CMB power spectrum, for e.g. by the use of spline-fitting, where the smoothed power 
spectrum has been used for visual aids (Oh et al. 1999)! 

In this paper the sparsity of the CMB power spectrum is used as a key ingredient in order to estimate the theoretical 
power spectrum without having to know the cosmological parameters; this estimate will not belong to a set of possible 
theoretical power spectra (i.e. all C[£] that can be obtained by CAMB 2 by varying the cosmological parameters). Instead, 
such an estimation should be useful for other applications, such as: 

• Monte Carlo: we may want to make Monte Carlo simulations in some applications without assuming the cosmological 
parameters. 

• Wiener filtering: Wiener filtering is often used to filter the CMB map and it requires the theoretical power spectrum 
as an input. We may not want to assume any cosmology at this stage of the processing. 

• Some estimators (weak lensing, ISW, etc.) require the theoretical power spectrum to be known. Using a data-based 
estimation of the theoretical C[£] could be an interesting alternative, or at least a good first guess in an iterative 
scheme where the theoretical C[£) is required to determine the cosmological parameters. 

Paper content 

Section 2 introduces the concepts of sparse representation and its applications to the CMB power spectrum. In Section 3 
we explain how sparsity is used to propose an estimator of the theoretical power spectrum. Experimental results are 
described and discussed in Section 4, and Section 5 is devoted to a discussion and comparison to the moving average 
estimator. In Section 6, conclusions are drawn and potential perspectives are stated. 



2 CAMB solves the Boltzmann equations for a cosmological model set out by the given cosmological parameters. 
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2. Sparsity of the CMB Power Spectrum 

2.1. A brief tour of sparsity 

A signal X — (X[l], . . . ,X[N]) considered as a vector in R , is said to be sparse if most of its entries are equal to zero. If 
k number of the N samples are not equal to zero, where k <C N, then the signal is said to be fc-sparse. In the case where 
only a few of the entries have large values and the rest are zero or close to zero the signal is said to be weakly sparse 
(or compressible). With a slight abuse of terminology, in the sequel, we will call compressible signals sparse. Generally 
signals are not sparse in direct space, but can be sparsified by transforming them to another domain. For example, sin(x) 
is 1-sparse in the Fourier domain, while it is clearly not sparse in the original one. In the so-called sparsity synthesis 
model, a signal can be represented as the linear expansion 

T 

X = $a = Y J <t>ia[i\ » (7) 

where a[i] are the synthesis coefficients of X, $ = (<fii, . . . , 4>t) is the dictionary, and <pi are called the atoms (elementary 
waveforms) of the dictionary <I>. In the language of linear algebra, the dictionary $ is a N x T matrix whose columns are 

the atoms normalized, supposed here to be normalized to a unit ^ 2 - n orm, i.e. Vi £ H&Ha = En=i I'M 71 ] 1 2 = I 3 - A 

function can be decomposed in many dictionaries, but the best dictionary is the one with the sparsest (most economical) 
representation of the signal. In practice, it is convenient to use dictionaries with fast implicit transform (such as Fourier 
transform, wavelet transform, etc.) which allow us to directly obtain the coefficients and reconstruct the signal from 
these coefficients using fast algorithms running in linear or almost linear time (unlike matrix- vector multiplications). 
The Fourier, wavelet and discrete cosine transforms provide certainly the most well known dictionaries. A comprehensive 
account on sparsity and its applications can be found in the monograph (Starck ct al. 2010). 



2.2. Which Dictionary for the Theoretical CMB Power Spectrum? 

We investigate the sparsity of the CMB power spectrum in two different dictionaries, both having a fast implicit transform: 
the Wavelet Transform (WT) and the Discrete Cosine Transform (DCT). 

Figure 1 shows an angular power spectrum (calculated by CAMB (Lewis et al. 2000) with WMAP7 (Larson et al. 
2010) parameters) along with the DCT- and WT-reconstructed power spectra with a varying fraction of the largest 
transform coefficients retained in the reconstruction. The inner plots show the difference between the actual power 
spectrum and the reconstructed ones. It can be seen that with only a few percentage of the coefficients the shape of the 
power spectrum is correctly reconstructed in both dictionaries. The height and the position of the peaks and troughs 
are of great importance here as the estimation of the cosmological parameters heavily relies on these characteristics of 
the power spectrum. The best domain would be the one with the sparsest representation and yet the most accurate 
representation of the power spectrum. Let C[£]^ M ' be its best M-term approximation, i.e. obtained by reconstructing 
from the M-largest (in magnitude) coefficients of C[£] in a given domain. To compare the WT and DCT dictionaries, we 
plot the resulting non- linear approximation (NLA) error curve in Figure 2, which shows the reconstruction error Em as 
a function of M, the number of retained coefficients; 



Em — 



C[t]-C[lf M) 



I CM || 



x 100 . (8) 



As M increases we get closer to the complete reconstruction and the error reaches when all the coefficients have been 
used. Usually the domain with the steepest E M curve is the sparsest domain. In this case though both dictionaries 
have very similar behaviors. There is only a small window in the coefficients for which DCT does a better job than 
WT. However, DCT flattens after using ~ 1% of the coefficients and does not improve the reconstruction until a big 
proportion of the coefficients have been used. 

Both dictionaries seem to suffer from boundary issues at low and high £s. This can be solved for high is as one 
can always perform the reconstruction beyond the desired I. For low £s it can be solved by different means, such as 
extrapolation of the spectrum. Note that the boundary issues are more severe in the DCT domain than WT; this is due 
to the fact that DCT atoms are not compactly supported. 

Next we investigate the sparsity of a set of realized spectra in the two dictionaries. We simulate 100 maps from the 
theoretical power spectrum used previously and estimate their power spectra using equation 4 . As before, we decompose 
each realization in the DCT and WT dictionaries and reconstruct keeping increasing fractions of the largest coefficients. 
At this stage, it is important to note that, as we are dealing with the empirical power spectrum, we are no longer in 
an approximation setting but rather in an estimation one. Indeed, the empirical power spectrum can be seen as a noisy 
version of the true one. Intuitively, reconstructing from a very small fraction of high coefficients will reject most of the 
noise (low estimator variance) but at the price of retaining only a small fraction of the true spectrum coefficients (large 
bias). The converse is true when a large proportion of coefficients is kept in the reconstruction. Therefore, there will exist 



The !p-norm of a vector X, p > 1, is defined as \\X\\ = (53, l-^M| J> ) 1 ■> with the usual adaptation H-Xll^ = max;X[i] 
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Fig. 1. A theoretical CMB power spectrum along with the reconstructed power spectra, using the DCT and WT dic- 
tionaries. The panels show the reconstructions for different fractions of the coefficients used. The inner plots show the 
differences between the actual and the reconstructed power spectra. Both dictionaries suffer from boundary effects, but 
this is more severe for DCT as the corresponding atoms are not compactly supported. It is worth mentioning that the 
power spectrum that is decomposed onto the two dictionaries is in the form £(£ + l)C{£]/2ir. 



a threshold value that will entail a bias-variance tradeoff, hence minimizing the estimation risk. This is exactly the idea 
underlying thresholding estimators in sparsifying domains. 

This discussion is clearly illustrated by the inner plots of Figure 3, which shows the normalized mean-square error 
(NMSE) defined as 

c[tf M) 



NMSE 



CI 



ran 



x 100 



(9) 



as a function of the fraction of coefficients used in the reconstruction. The error is large when only a few coefficients are 
used. As more coefficients are included, one starts to recover the main (i.e. the general shape of the spectrum) features of 
the power spectrum. With more coefficients, more noise enters the estimation and the error increases again. The NMSE 
curve shows a clear minimum at which the underlying true power spectrum is best recovered. 

Despite the differences in the performance of the two dictionaries, the minima of the NMSE are around the same 
proportions of the coefficients. This is because the NMSE reflects a global behavior. On the one hand, although the DCT 
can recover the features of the spectrum correctly, it is less smooth than WT. Conversely, the WT cannot reconstruct 
the proper shape of the power spectrum, but provides a smoother estimate. 

Figures 4 and 5 show the same results for an average over the 100 realizations. It can be seen that on average DCT 
does a great job in the recovery of the features of the CMB spectrum. Indeed the minimum of NMSE curve is at a lower 
proportion of coefficients for DCT than WT. This is because the small noisy features of the DCT-reconstructed spectra 
cancels out in the averaging, while the wrong recovery of the shape of the spectrum by WT does not. In a nutshell, DCT 
seems to do a better job in reconstructing the true underlying CMB power spectrum from its realizations. 

To summarize, from the above discussion, we conclude the following: 
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Fig. 2. Non-Linear Approximation (NLA) error curves for the two dictionaries. Below 1% the DCT curve is dropping 
faster, which means it is doing a better job. However, past ~ 2% the DCT curve flattens off while WT decreases to ~ 
very quickly. 

• the CMB power spectrum is very sparse in both the DCT and WT dictionaries, although their sparsifying capabilities 
are different; 

• DCT recovers global features of spectrum (i.e. the peaks and troughs) while WT recovers localized features; 

• in the case of realizations, WT recovers more localized (noisy) features than the global ones, while the DCT concen- 
trates on the global features. 

In the next section, these complementary capabilities of the DCT and WT transforms will be combined to propose a 
versatile way for adaptively estimating the theoretical power spectrum from a single realization of it. 

3. Sparse Reconstruction of the Theoretical Power Spectrum 

Let's start with the simple model where the observed signal Y is contaminated by a zero-mean white Gaussian noise, 
Y = X + e, where X is the signal of interest and e ~ Af(0,a 2 ). Sparse recovery with an analysis-type sparsity prior 
amounts to finding the solution of the following problem: 

mm^X^ s.t. \\Y-X\\ 2 <5, (10) 

where <fr T X represent the transform coefficients of X in the dictionary $, and 5 controls the fidelity to the data and 
obviously depends on the noise standard deviation a. 

Let's now turn to denoising the power spectrum from one empirical realization of it. In this case, however, the noise 
is highly non-Gaussian and needs to be treated differently. Indeed, as we will see in the next section, the empirical power 
spectrum will entail a multiplicative x 2 -distributed noise with a number of degrees of freedom that depends on I. That 
is, the noise has a variance profile that dependents both on the true spectrum and t. We therefore need to stabilize the 
noise on the empirical power spectrum prior to estimation, using a Variance Stabilization Transform (VST). Hopefully, 
the latter will yield stabilized samples that have (asymptotically) constant variance, say 1, irrespective of the value of 
the input noise level. 



(i 
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Fig. 3. A simulated CMB power spectrum along with the reconstructed spectra, using the DCT and WT dictionaries. 
The black solid line is the true underlying power spectrum from which the simulations were made. The blue and red dots 
show the simulated and the reconstructed power spectra respectively. With only 1% of the coefficients, DCT can recover 
the input power spectrum (i.e. the black solid line) very well, recovering the peaks and troughs accurately. Unlike DCT, 
WT seems to have difficulties in recovering the peaks and troughs. The inner plots shows the NMSE curves. 

3.1. Variance Stabilizing Transform 

In the statistical literature the problem of removing the noise from an empirical power spectrum goes by the name of 
periodogram denoising (Donoho 1993). In (Komm et al. 1999), approximating the noise with a correlated Gaussian noise 
model, a threshold was derived at each wavelet scale using the MAD (Median of Absolute Deviation) estimator. A more 
elegant approach was proposed in (Donoho 1993; Moulin 1994), where the so-called Wahba VST was used. This VST is 
dehned as: 

T(V) = (logV + 7 )— , (11) 

where 7 = 0.57721... is the Euler-Mascheroni constant. After the VST, the stabilized samples can be treated as if the 
noise contaminating them were white Gaussian noise with unit variance. 

We will take a similar path here, generalizing the above approach to the case of the angular power spectrum. Indeed, 
from (4), one can show that under mild regularity assumptions on the true power spectrum, 

C[£] A C[£]Z[£], where W > 2, 2LZ[£] ~ X Il> L = 2£ + l. (12) 

-4- means convergence in distribution. From (12), it is appealing then to take the logarithm so as to transform the 
multiplicative noise Z into an additive one. The resulting log-stabilized empirical power spectrum reads 

C s [£] := T e (C[£}) = \ogC[£] - p L = logC[£] + r,{£] . (13) 

where r)[£] := logZ[£] — L = 2£ + 1. Using the asymptotic results from (Bartlett & Kendall 1946) on the moments 
of log— x variables, it can be shown that \x L — ij) (L) — logL, E(r?[£]) = and a\ — Var [q[£]] = ipi (L), where ip m (t) is 

the standard polygamma function, ip m (t) = £ tm +i logr(i). 
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Fig. 4. The average of the DCT-reconstructed power spectra when different fraction of most significant coefficients have 
been used in the reconstruction. The inner plot shows the NMSE curve for this average. The minimum of the curve is 
at less than 10% of the coefficients, meaning that the true CMB spectrum can be recovered with less than 10% of the 
coefficients while ensuring a good bias-variance tradeoff. 



We can now consider the stabilized C s [£] as noisy versions of the log(7[£], where the noise is zero-mean additive and 
independent. Owing to the Central Limit Theorem, the noise tends to Gaussian with variance <j\ as £ increases. At low 
£, normality is only an approximation. In fact, it can be show that the noise rj[£] has a probability density function of 
the form 

p,W = J^ex P [L(^ + ML - e ^)] , (14) 

which might be used to estimate the thresholds in the wavelet domain. 

In order to standardize the noise, the VST (13) will be slightly modified to the normalized form 

C s [£] := T e {C[£}) = logC ^~^ L = X s [£] + e\l] . (15) 

where now the noise e[£] is zero-mean (asymptotically) Gaussian with unit variance, and X s [£] := log C[£]/<tl- It can be 
checked that the Wahba VST (11) is a specialization of (15) to L = 0. 

In the following, we will use the operator notation T(X) for the VST that applies (15) entry-wise to each X[£], and 
H{X) its inverse operator, i.e. U(X) := (H e (X[£\)) t with TZ e {X[£}) = cxp{a L X[£}). 



3.2. Signal detection in the wavelet domain 

Without of loss of generality, we restrict our description here to the wavelet transform. The same approach applies to 
other sparsifying transforms, e.g. DCT, just as well. 

In order to estimate the true CMB power spectrum from the wavelet transform, it is important to detect the wavelet 
coefficients which are "significant" , i.e. the wavelet coefficients which have an absolute value too large to be due to noise 
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Fig. 5. Same as Figure 4, but for the WT domain. It seems that WT cannot do a great job for simulated power spectra, 
compared to the DCT domain; the minimum of the NMSE curve is at more than 35% — while the curve is pretty much 
flat after ~ 20% anyway. 



(cosmic variance + instrumental noise). Let Wj[£] the wavelet coefficient of a signal Y at scale j and location I. We define 
the multiresolution support M of Y as: 

_ J 1 if W j [£} is significant, 
^ J ~\0 if Wj [£\ otherwise. [ ' 

For Gaussian noise, it is easy to derive an estimation of the noise standard deviation <jj at scale j from the noise standard 
deviation, which can be evaluated with good accuracy in an automated way (Starck & Murtagh 1998). To detect the 
significant wavelet coefficients, it suffices to compare the wavelet coefficients in magnitude to a threshold level tj. 

This threshold is generally taken to be equal to Kdj, where k ranges from 3 to 5. This means that a small magnitude 
compared to the threshold implies that the coefficients is very likely to be due to noise and hence insignificant. Such a 
decision rule corresponds to the hard-thresholding operator 

if | re j [€] | > tj then Wj[£] is significant , 

if < tj then Wj [£] is not significant. ^ ' 

To summarize, The multiresolution support is obtained from the signal Y by computing the forward transform 
coefficients, applying hard thresholding, and recording the coordinates of the retained coefficients. 



3.3. Power Spectrum Recovery Algorithm 

Let's now turn to the adaptive estimator of the true CMB power spectrum C[£] from its empirical estimate C[£}. As we 
benefit from the (asymptotic) normality of the noise in the stabilized samples C s [£] in (15), we are in position to easily 
construct the multiresolution support M of C s as described in the previous section. Once the support M of significant 
coefficients has been determined, our goal is reconstruct an estimate X of the true power spectrum, known to be sparsely 
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represented in some dictionary $(regularization), such that the significant transform coefficients of its stabilized version 
reproduce those of C s (fidelity to data). Furthermore, as a power spectrum is a positive, a positivity constraint must 
be imposed. These requirements can be cast as seeking an estimate that solves the following constrained optimization 
problem: 

T 1 "* 7 *" 1 **' {mI%tT(X))=Mq(*t C *) ' (18) 

where stands for the Hadamard product (i.e. entry-wise multiplication) of two vectors. This problem has a global 
minimizer which is bounded. However, beside non-smoothness of the Zi-norm and the constraints, the problem is also 
non-convex because of the VST operator T. It is therefore far from obvious to solve. 

In this paper we propose the following scheme which starts with an initial guess of the power spectrum X^ = 0, and 
then iterates for n — to ]V max — 1, 

x = n(r(x {n) ^ +$m© ($ T (c s -t 

x (n+l) = v+ ^ ST Xn (^ T X)) , 

where V+ denotes the projection on the positive orthant and guarantees non-negativity of the spectrum estimator, 
STa„(u>) = (ST\ n (w[i])) i is the soft-thresholding with threshold A n that applies term-by-term the shrinkage rule 

ST (w\t\) = { sign ( w H)(l w HI ~ A «) if A ™ ' (20) 

X " { [l> \0 otherwise. 1 ' 

Here, we have chosen a decreasing threshold with the iteration number n, X n = (N max — n)/{N max — 1). More details 
pertaining to this algorithm can be found in Starck et al. (2010) 

Algorithm 1 summarizes the main steps of the sparse denoising algorithm. A similar approach was proposed in (Starck 
et al. 2009; Schmitt et al. 2010) for Poisson noise removal in 2D and 3D data sets. 



Algorithm 1: TOUSI Power Spectrum Smoothing 



Require: 

Empirical power spectrum C, 



9 

10 



Threshold k (default value is 5). 
Detection 

Compute C s using (15). 

Compute the decomposition coefficients W of C s in $, W = $ T C". 

Compute the support M from W with the threshoid k, assuming standard additive white Gaussian noise. 

Estimation 

Initialize X m = 0. 

for n = to 7V max — 1 do 

x = ii (t +M0 ($ T (c s - r (x (n) )))). 

\ _ JV max -(n + l) 
n + 1 _ JVmax-1 ' 

end for 

Return: The estimate X = X (JVmax) . 



3.4. Instrumental Noise 

In practice, the data are generally contaminated by an instrumental noise, and estimating the true CMB power spectrum 
C[£] from the empirical power spectrum C[£] requires to remove this instrumental noise. The instrumental noise is 
assumed stationary and independent from the CMB. We will also suppose that we have access to the power spectrum 
of the noise, or we can compute the empirical power spectrum Sn[£] of at least one realization, either from a JackKnife 
data map or from realistic instrumental noise simulations. The above algorithm can be adapted to handle this case after 
rewriting the optimizing problem as follows: 



(21) 
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Thus, (19) becomes 

X = k(t (x™ + SjA + $M ($ T (c s - T (xW + S N y\ )) - S N 
X (n+i) = v+ U ST An ($ T X)) . 

Algorithm 1 can be modified accordingly. 

3.5. Combining Several Dictionnaries 

We have seen in Section 2.2 that the WT and DCT dictionaries had complementary benefits. Indeed each dictionary 
is able to capture well features with shapes similar to its atoms. More generally, assume that we have D dictionaries 
, <f>£). Given a candidate signal Y, we can derive a support Md associated to each dictionary $d, for d £ {1, • • ■ , D}. 
The optimization problem to solve now reads 

Again, this is a challenging optimization problem. We propose to attack it by applying successively and alternatively 
(22) on each dictionary $d- Algorithm 2 describes in detail the different steps. 



Algorithm 2: TOUSI Power Spectrum Smoothing with D dictionaries 



Require: 

Empirical power spectrum C, D dictionaries $i, $.d, noise power spectrum Sn, 
Number of iterations iV max , 
Threshold k (default value is 5). 
Detection 



1 

2 
3 

1 
5 
6 
7 
8 

9 

10 
11 
12 
13 
11 



Compute C s using (15). 

For all d, compute the decomposition coefficients Wd of C s in $cj, Wd = $JC S . 

For all d, compute the support Md from Wd with the threshold k, assuming standard additive white Gaussian noise. 

Estimation 

Initialize X (0) = 0, 

for n = to -/V max — 1 do 

Z d =X&>. 

for d — 1 to D do 



Z = TZ (T (Z d + S N ) + $ d M ($ d T (<7 S - T (z d + &v)))) - S N - 
Z d +i=V+ ($ d ST A „ ($/£)). 



end for 

V-(«+l) _ ^-0 + ! 

\ _ JVmax-(n + l) 

An+1 _ JV max -l ' 

end for 

Get the estimate X = X (N " 



4. Application: Monte Carlo Simulations 

We simulate 100 maps from a theoretical CMB power spectrum that is calculated by CAMB. The power spectra of these 
maps are equivalent to 100 realizations of the true CMB power spectrum — This realized spectrum is what we have 
access to in reality. Each of these 100 simulated spectra are run through the TOUSI algorithm, with the aim of recovering 
the theoretical spectrum from which these 100 spectra were simulated (i.e. the one that was calculated by CAMB). 

Figure 6 shows an example of this; the empirical power spectrum of one realization (blue dots) that was fed into the 
TOUSI algorithm, the average of the 100 estimated power spectral using TOUSI (red line) and the input theoretical 
spectra (black line). The black line is the input theoretical spectrum that was calculated by CAMB, which is what we 
are trying to recover. The reconstruction of the peaks and troughs of the power spectrum by this algorithm is very 
impressive. This is very important as these features define the cosmological parameters. To further check the accuracy of 
these reconstructed spectra, we estimate a set of cosmological parameters from these spectra, using CosmoMC (Lewis & 
Bridle 2002). First, a set of cosmological parameters are estimated from the 100 simulated spectra. The results are shown 
in Figure 4 as black solid lines. Then a 'mean' reconstructed spectrum is calculated by averaging the 100 reconstructed 
spectra. This, in principle (i.e. if the algorithm has worked), should be an estimation of the true input spectrum with 
the same characteristics and the same cosmological parameters. To test this, we use this 'mean' reconstructed spectrum 
to simulate another 100 maps and then 100 spectra. These simulated spectra are run through CosmoMC to estimate the 
same set of cosmological parameters. These are shown as red lines in Figure 4. It can be seen that TOUSI algorithm can 
reconstruct the true underlying power spectrum with great accuracy in the cosmological parameters. 
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Fig. 6. The theoretical CMB power spectrum (black line), the empirical power spectrum of one realization (blue dots) 
and the avegared of estimated power spectra (red line) using TOUSI algorithm. The inner plots show a zoomed-in version. 

4.1. Data with Instrumental Noise 

Here we present the performance of the TOUSI algorithm in the presence of instrumental noise. The noise maps were 
simulated using a theoretical (PLANCK level) noise power spectrum. They were added to the CMB maps simulated 
previously and the power spectra of the combined maps were estimated using equation 4. 

Figure 8 shows the reconstruction of the theoretical CMB spectrum in the presence of noise. The blue dots show the 
empirical power spectrum of one realization having instrumental noise. Yellow dots show the estimated power spectrum 
of one of the simulated noise maps. Green dots show the the spectrum with the noise power spectrum removed. The 
black and red solid lines are the input and reconstructed power spectra respectively. The theoretical power spectrum can 
be reconstructed up to the point where the structure of the power spectrum has not been destroyed by the instrumental 
noise. In our case, having PLANCK level noise, this goes to I up to 2500. It can be seen that TOUSI can do a great job 
in reconstructing the input power spectrum even in the presence of instrumental noise. 

4.2. Test on WMAP7 power spectrum 

We test our algorithm on real data. Figure 9 shows the application of our technique to the WMAP7 power spectrum. 
The method works really well up to I of ~ 800. After this point the instrumental noise becomes so dominant that the 
features of the spectrum is washed out and cannot be recovered. 



5. Sparsity versus Averaging 

A very common approach to reduce the noise on the power spectrum is the moving average filter, i.e. average values in 
a given window, 



C A [b] = 



1 

b(b+ l)oj b 



(24) 



i=b- 
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Fig. 7. Cosmological parameters estimated from the true CMB power spectrum (black line) and the mean of reconstructed 
power spectra (red line). The dashed line is the true input parameters, i.e. the ones used to calculate the theoretical 
power spectrum using CAMB. 

and the window size w& is increasing with i. Here, we use window sizes of {1, 2, 5, 10, 20, 50, 100} respectively for I ranging 
from {2, 11, 31, 151, 421, 1201, 2501} to {10, 30, 150, 420, 1200, 2500, 3200}, which have also been used in the framework of 
the PLANCK project in (Leach ct al. 2008). 

Figure 10 shows a reconstruction of the power spectrum for TOUSI versus averaging. From the NMSE curve it is 
clear that our algorithm is much more efficient. Figure 11 shows the average error the 100 realizations as a function of I 

too 

E ® = looEllCM-CiMlh, (25) 

i—l 

where Ci is the estimated power spectrum from the z-th realization. We display the errors for the spectra estimated by 
the empirical estimator (the realization, black dotted line), the averaging estimator (red dashed line) and TOUSI (solid 
blue line). The cosmic variance is over-plotted as a solid black line. We can see that the expected error is highly reduced 
when using the sparsity-based estimator. 

6. Conclusion 

Measurements of the CMB anisotropics are powerful cosmological probes. In the currently favored cosmological model, 
with the nearly Gaussian-distributed curvature perturbations, almost all the statistical information are contained in 
the CMB angular power spectrum. In this paper we have investigated the sparsity of the CMB power spectrum in two 
dictionaries; DCT and WT. In both dictionaries the CMB power spectrum can be recovered with only a few percentages 
of the coefficients, meaning the spectrum is very sparse. The two dictionaries have different characteristics and can 
accommodate reconstructing different features of the spectra; The DCT can help recover the global features of the 
spectrum, while WT helps recover small localized features. The sparsity of the CMB spectrum in these two domains 
has helped us develop an algorithm, TOUSI, that estimates the true underlying power spectrum from a given realized 
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Fig. 8. Power spectrum estimation in the presence of instrumental noise. The blue dots show the empirical power 
spectrum of one realization having instrumental noise. Yellow dots show the estimated power spectrum of one of the 
simulated noise maps. Green dots show the the spectrum with the noise power spectrum removed. The black and red 
solid lines are the input and reconstructed power spectra respectively. The inner plots show a zoomed-in version. 



spectrum. This algorithm uses the sparsity of the CMB power spectrum in both WT and DCT domains and takes the 
best from both worlds to get a highly accurate estimate from a single realization of the CMB power spectrum. This 
could be a replacement for CAMB in cases where knowing the cosmological parameters is not necessary. The developed 
IDL code will be released with the next version of ISAP (Interactive Sparse astronomical data Analysis Packages) via 
the web site: 



http : // j starck . free . f r/isap . html 
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Fig. 10. Mean denoised spectra with Wavelets (blue solid line) and averaging (red dashed line) from 100 realizations. 
The inner plot shows the normalized error for both dictionaries. 
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Fig. 11. Mean error for the 100 realizations, for the realizations (black dotted line), the averaging denoising (red dashed 
line) and the sparse wavelet filtering (blue solid line). The inner plot shows a zoom between / = 2000 and / = 3000. 



